function [ROC] = kneu_ROC(testdat,testlabels,W,b0,do_plot)
% [test_ROC train_ROC wdist] = kneu_ROC(testdat,testlabels,W,b0,do_plot)
% Example using the spider toolbox :
% d=data([rand(500,2);rand(500,2)+0.8],[-ones(500,1); ones(500,1)]);
% dt=data([rand(500,2)+0.25;rand(500,2)+0.55],[-ones(500,1); ones(500,1)]');
% alg=svm({'C=1'});
% [tr,alg]=train(alg,d);
% %visualize data:
% w=get_w(alg);
% plot(d.X(d.Y<0,1),d.X(d.Y<0,2),'b.',d.X(d.Y>0,1),d.X(d.Y>0,2),'r.',dt.X(dt.Y<0,1),dt.X(dt.Y<0,2),'co',dt.X(dt.Y>0,1),dt.X(dt.Y>0,2),'mo',[0 2],[-alg.b0 -2*w(1)-alg.b0]/w(2),'k-');
% ROC = kneu_ROC(dt.X,dt.Y,get_w(alg),alg.b0,true);


if nargin<5,
    do_plot=true;
end

test_size = size(testdat,1); % number of train samples

x_dot_w_test = testdat * (W'*ones(1,test_size));
x_dot_w_test = (x_dot_w_test .* eye(test_size)) * ones(test_size,1);
% get separate indices
idx1 = find(testlabels>0);
idx2 = find(testlabels<0);
N1 = length(idx1); % maximum of possible true elements
N2 = length(idx2); % maximum of possible false elements
% distances to hyperplane (sorted)
[wdist sortIdx] = sort((x_dot_w_test+b0)/norm(W));
sortedTrueLabels = testlabels(sortIdx)';
wdist1_test = sort((x_dot_w_test(idx1)+b0)/norm(W));
wdist2_test = sort((x_dot_w_test(idx2)+b0)/norm(W));

% ROC 
ROC.TP = zeros(1,test_size); % true positives
ROC.FP = zeros(1,test_size); % false positives
ROC.FN = zeros(1,test_size); % false negatives
ROC.TN = zeros(1,test_size); % true negatives
ROC.b = zeros(1,test_size); % shifted criterion
for k=1:test_size,
    ROC.TP(k)=sum(sortedTrueLabels>0 & wdist>=wdist(k));
    ROC.FP(k)=sum(sortedTrueLabels<0 & wdist>=wdist(k));
    ROC.FN(k)=sum(sortedTrueLabels>0 & wdist<wdist(k));
    ROC.TN(k)=sum(sortedTrueLabels<0 & wdist<wdist(k));
    ROC.b(k)=b0+wdist(k)*norm(W);
end
ROC.recall=[ROC.TP./(ROC.TP+ROC.FN); ROC.TN./(ROC.TN+ROC.FP)];
ROC.precision=[ROC.TP./(ROC.TP+ROC.FP); ROC.TN./(ROC.TN+ROC.FN)];
ROC.accuracy=(ROC.TP+ROC.TN)./test_size;
ROC.TP=ROC.TP/N1;
ROC.FP=ROC.FP/N2;
ROC.FN=ROC.FN/N1;
ROC.TN=ROC.TN/N2;
ROC.AUC = sum(ROC.TP.*diff([0 1-ROC.FP]));
[ROC.opt_b,optPos]=max(abs(ROC.accuracy));

if do_plot,    
    % plot results
    figure;    
    %distribution along hyperplane normal
    nbins = max(10,sqrt(size(testdat,1)));
    hold on
    [N,X]=hist(wdist1_test,nbins);
    bH=bar(X,N,'b');
    set(get(bH,'children'),'facealpha',0.5);
    [N,X]=hist(wdist2_test,nbins);
    bH=bar(X,N,'r');
    set(get(bH,'children'),'facealpha',0.5);
    xlabel('Distance to Hyperplane')
    ylabel('Number of Samples')
    lH=legend('Cond1', 'Cond2', 'Location','NorthEastOutside');
    box off;
    legend boxoff;
    set(lH,'fontsize',20)    
    hold off
    % ROC train: call function with train data as test data    
    % ROC test
    figure;
    hold on;
    plot(ROC.FP,ROC.TP,'linewidth',2);
    ylim([0 1.01]);
    [dummy,minPos]=min(abs(ROC.b-b0));
    plot(ROC.FP(minPos),ROC.TP(minPos),'go') % actual criterion
    plot(ROC.FP(optPos),ROC.TP(optPos),'ro') % optimal criterion
    title(sprintf('Test data; Area under ROC: %1.4f',ROC.AUC));
    xlabel('False Alarms');
    ylabel('Hits');
    hold off
end
